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SUMMARY 

Some recent developments in acceleration of convergence methods for vector 
sequences are reviewed. The methods considered are the minimal polynomial 
extrapolation, the reduced rank extrapolation, and the modified minimal poly- 
nomial extrapolation. The vector sequences to be accelerated are those that 
are obtained from the iterative solution of linear or nonlinear systems of 
equations. The convergence and stability properties of these methods as well 
as different ways of numerical implementation are discussed in detail. Based 
on the convergence and stability results, strategies that are useful in 
practical applications are suggested. Two appl i cations to computational fluid 
mechanics involving the three-dimensional Euler equations for ducted and 
external flows are considered. The numerical results demonstrate the useful- 
ness of the methods in accelerating the convergence of the time-marching 
techniques in the solution of steady state problems. 


1. INTRODUCTION 

Let R^ denote the d-dimensional Cartesian space and let Q be a subset 
of R^. Denote the boundary of Q by 3Q. Let u(t) be the solution to the 
initial boundary value problem (IBVP) 



du 

dt + 

4>(u) = 

0 in Q, 

A(u) 

= 0 on 

3Q 

(boundary condition). 

M(u) = 

0 at 

ct- 

II 

O 

(initial condition). 


Here <|>, A, and M are linear or nonlinear operators (differential, integral, 
etc.), and t denotes time. 
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Assume that the IBVP in equation (1.1) has a time-independent (steady 
state) solution that we shall denote u*. Then u* is the solution to the 
boundary value problem (BVP) 


<(>( u* ) =0 in Q, 

A(u* ) =0 on 8Q (boundary condition). 


( 1 . 2 ) 


Since u* = 1 i m u(t), one of the widely used techniques for determining 

t">co 

u* has been one in which the IBVP in equation (1.1) is solved by a time- 
marching technique. Specifically, equation (1.1) is discretized to give the 
iterative scheme 


u n+1 = i|;( u n ) , n = 0,1, . . ., (1.3) 

where u n is the (discrete) approximation to u(nAt) and At > 0 is the time 

increment. The initial approximation u® needs to be provided in an 
appropriate way. Obviously, when the discretization scheme in equation (1.3) 
is stable lim u n = u, where u is the fixed point of the function <|r(z) and 

n-»oo 

in addition is a discrete approximation to u*. In practice, when 
| | u n+ l - u n | | / | | u 1 - u^| | < e, for some prescribed e > 0, u n is taken to be 

an acceptable approximation to u. Note that u n+ ^ - u n = i|/(u n ) - u n , which 

is the residual of <|>(z) - z for z = u n , and ||«|| denotes some appropriate 
vector norm. 

When the discretization scheme in equation (1.3) is stable, for n suffi- 
ciently large, u n is close to u, hence 

U n+ 1 = <|»(u n ) = >|/(u) + f'(u)(u n - U) + e n = Au n + b + e n , (1.4) 

where A = y'(u) is the Jacobian of f at u, b = y(u) - <j»‘(u)0, and 
e n satisfies ||e n || < C||u n - u | | 2 for some fixed constant C. Obviously 
A and b are independent of n. That is to say, the u n , for sufficiently 
large n, satisfy (approximately) 

u n+ l =; Au n + b. (1.5) 


Thus 


1 1 u n — u 1 1 = 0{p n Ct'(u)]} as n -> ®, (1.6) 

where, for any matrix B, p(B) denotes its spectral radius. 

Normally p[f(u)] is very close to 1, although it is_strictly less than 
1. In addition, in some cases it can be shown that p[\|»'(u)] tends to 1 as the 
meshsize of the discretization tends to zero. This causes the sequence u^, u 1 , 
u^, .... to converge to u very slowly. One simple way to overcome this 
problem is by use of convergence acceleration (or equivalently extrapolation) 
methods that do not require that changes be made in the basic iterative scheme 
of equation (1.3). In the next section we shall describe three such methods, 
namely, the minimal polynomial extrapolation (MPE) of Cabay and Jackson 
(ref. 3), the reduced rank extrapolation (RRE) of Eddy (ref. 5) and Mesina 
(ref. 12), and the modified minimal polynomial extrapolation (MMPE) of Sidi, 
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Ford, and Smith (ref. 19). A slightly different version of RRE was earlier 
proposed by Kaniel and Stein (ref. 11). All three methods operate on the given 
vector sequence u^, , u 2 , . . . , only, irrespective of how this sequence 

is obtained (or, in the context of the discretization scheme (eq. (1.3)), 
irrespective of what y is). This is an important property of these methods 
as it allows them to be applied in conjunction with iterative methods for both 
linear and nonlinear problems with the same ease. We must add, however, that 
the rates of acceleration that can be achieved by using these methods depend 
on the sequence in consideration. Hence, in the context of equation (1.3), it 
is the structure of t that determines the rates of acceleration. 

In the next section we shall give a brief description of MPE, RRE, and 
MMPE. These descriptions are based on the developments of the survey paper of 

Smith, Ford, and Sidi (ref. 20) and of ref. 19 and of Sidi (ref. 15). Further 

generalizations that are proposed in reference 19 will also be mentioned. The 
convergence and stability properties of these methods have been analyzed in 
references 15 and 19, Sidi and Bridger (ref. 18), and Sidi (ref. 16). Some of 

these results will be reviewed in section 3. In section 4 we shall present 

some applications to certain three-dimensional fluid mechanics problems. 

Extrapolation methods have recently been used by several authors in 
computational fluid mechanics applications, see, for example, Hafez et al . 

(ref. 7), Wong and Hafez (ref. 22), Wigton, Yu, and Young (ref. 21), Jespersen 
and Buning (ref. 10), and Reddy and Jacocks (ref. 14). One of the methods 
described in the present work, namely, MPE, with some variation, is employed 
in references 10 and 14. 

Before closing this section we mention that the survey paper (ref. 20) 
discusses in detail MPE and RRE, as well as the well known scalar epsilon 
algorithm (SEA) of Wynn (ref. 23) and its two vector versions, the vector 
epsilon algorithm (VEA) of Wynn (ref. 24) and the topological epsilon algorithm 
(TEA) of Brezinski (ref. 2). SEA and VEA have been used in accelerating the 
convergence of some numerical schemes in computational fluid mechanics, see, 
e.g., Hafez et al . (ref. 7). As explained in reference 20 and in Ford and Sidi 
(ref. 6), the epsilon algorithms are expensive both storagewise and timewise. 
They need about twice as much storage as MPE, RRE, or MMPE, and their vector 
operation counts are several times those of MPE, RRE, or MMPE. 


2. DESCRIPTION OF VECTOR EXTRAPOLATION METHODS 

We shall now give a brief description of MPE, RRE, and MMPE based on the 
developments in references 20, 15, and 19. 

Let xq , xi , X£, . . • , be a vector sequence in the complex N-dimensional 
Euclidean space C N . Assume that this sequence converges, and denote its limit 
by s. Using the vectors Xj only, MPE, RRE, and MMPE produce approximations 
to s, which are determined as follows: 

Def i ne 

u i = Ax} = Xi + i - Xi, w, = Au j = A 2 Xi, i = 0,1,2, ... . (2.1) 
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MPE - Minimal Polynomial Extrapolation 
Pick a positive integer k < N and form the matrix U by 


, u n ! u 


n+1 


! U n+k-l. 


( 2 . 2 ) 


and solve the overdetermined (and in general inconsistent) system of equations 


Uc = -u 


n+k’ 


(2.3) 


where c = (c Q , c 1 , . . . , c^)^ by linear least squares. With c Q , c 1 , 
. . . , c^ 1 determined, set c k = 1 , and let 


Y 
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j = 0,1, . . . , k. 


Finally set 


(2.4) 


; n k = ^ Y i 
n,K j=0 3 


n+j ’ 


where s , is the desired approximation to s. Note that 

D , 


£ Y-i = 1 

j-o 3 


(2.5) 


( 2 . 6 ) 


as is implied by equation (2.4), i.e., s nj k is a weighted "average" of the 
vectors x n , x n+ i , ...» x n+ k, in which’the weights are not necessarily real 
and nonnegative. 

If we define the inner product associated with C N to be 


where y = (V . • 
in C^, then the 


N T i 

(y,z) = y*z = yV , (2.7) 

i = l 

, y N ) and z = (z 1 , . . . , z N ) are arbitrary vectors 
are equivalently the solution to the normal equations 


k-1 

£q (u n+i’ U n+j )c j = _(u n+i’ u n+k )l 


i = 0,1 , 


k-1. (2.8) 


Finally, if we let U + be the generalized inverse of the matrix U, then 
the vector c is also given by 


c = -U + u 


n+k' 


(2.9) 
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RRE - Reduced Rank Extrapolation 


Pick a positive integer k < N and form the matrix W by 


H 



i i 

1 w 1 
: n + i ; 


i w n+k-l 


- 


( 2 . 10 ) 


and solve the overdetermined (and in general inconsistent) system of equations 

Wq = -u n , (2.11) 

where q = (q Q , q ] , . . . , q k _.,) T , b Y linear least squares. With q Q , q ] , 

. . . , q [( ._ 1 determined, set 


where s n>k 
The q. 


k-1 

S n,k ■ X n * £ «jVr 

is the desired approximation to s. 
are equivalently the solution to the normal equations 


( 2 . 12 ) 


k-1 

£ 0 <Vr V3X3 ■ - ( Vi' V- 1 ■ °-' 


k-1. (2.13) 


Also, if W + is the generalized inverse of the matrix W, then the vector q 
is given by 

q = -W + u n • (2.14) 

Interestingly enough, s n>k for RRE, just like s nj k f° r Mp E> can be 
expressed as in equation (2.5)! with equation (2.6) holding true. Specifi- 
cally, the qi and the corresponding yj for RRE are related by 

y 0 - ' - q o ; -"j-i - V \-\-r <2 - 15> 


MMPE - Modified Minimal Polynomial Extrapolation 

Pick a positive integer k < N and form the matrix U as in equa- 
tion (2.2) and "solve" the overdetermined system of equations in equation (2.3) 
for c = (cq, ci, ... , c k _])T, in the sense 

QUc = -Qufc, (2.16) 

where Q is a f i xed k x N matrix of full rank. Obviously the matrix QU of 
equation (2.16) is k x k. If qi, . . . , q k are the rows of the matrix Q, 
then equation (2.16) can also be expressed as 


k-1 

£ <V Vj )c 3 ’ - ( V u n+k> ’ 1 * 1 ± k - 


(2.17) 
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c.f. equation (2.8) for MPE. Now set = 1 and determine the y j , 
0 < j < k, by equation (2.4), and s n |<_, by equation (2.5). 


Determinantal representations for s n> k exist both for MPE, RRE, and 
MMPE, and they are of the form 


D( V x n+l 


• x n+k > 


*n,k 


D( 1 , 1 ..... 1 ) 


where D(<Jq, a 1 , . . . , a^) is the determinant 


(2.18) 


D(a 0> a 1 , . . . , o^) = 


J 0,0 


J 0,1 


0,k 


(2.19) 


and 


u 


k-1 ,0 


u 


k-1 ,1 


'k-1 


,k 




n+f 

for 

MPE, 


n+f 

for 

RRE, 

(2.20) 

n+f 

for 

MMPE. 



When the c. are vectors D(ct q , , .... a^) is also a vector, and is 
k 

defined to be £ a.N., where N. is the cofactor of a.. 
i=0 1 11 

The approximations s n ^ for all three methods are determined from the 
k + 2 vectors x n , x n+ i , , x n+ k, as can easily be seen from their 

definitions. 

As far as the implementation of MPE and RRE is concerned, a few alterna- 
tives exist: 

(1) Solution by least squares packages: Although feasible for small scale 

problems, this approach may become very costly for large scale problems in 
which N, the dimension of the vectors, is very large and k is not small. 

We may even run into storage problems, as the present least squares solvers 
require the matrix of the equations to be solved (U for MPE and W for RRE) 
to be stored in core memory. One way out of this limitation would be to modify 
these solvers so that this matrix can be stored in a secondary storage device 
like a disk, for example. Additional storage for the k + 2 vectors x n , 

. . . , x n+ k + i is also required. 
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(2) Solution of normal equations: Although the conditioning of this 

approach is less favorable than that for the previous one, the storage require- 
ments for it are much lower. Actually, with appropriate programming, we can 

see that only k + 2 vectors, namely x n , u n , u n+ i u n+ k, need to be 

saved. That this is so can be shown as follows: Define u^j E (u n+ i , u n+ j). 

(a) For MPE the system of equations in equation (2.8) is entirely given 
in terms of the U} y Once the c i , hence y^ have been obtained from equa- 
tions (2.8) and (2.4;, s n> k in equation (2.5) can be computed from 


k-1 

s n,k = *n + ^i u n+i ’ 


( 2 . 21 ) 


where the £< can be determined recursively from the y. through the 
relations J 


Sq ~ i Yq» ~ ~ Yj > j = • • • > k - 1, (2.22) 

or the relations 

e k-l = Y k : ^j-1 = Yj + Zy j = k - l,k - 2 0. (2.23) 


(b) For RRE the system of equations in equation (2.13) is determined solely 
in terms of the uyy since (w n+1 , w n+ j) = u i + i j + i + Oi j - Gi +1 j - uij + i- 
Once the qj have been determined from equation (2.13), s n ^ is obtained 
directly from equation (2.12). 


(3) Recursive computation: Both MPE and RRE can be implemented by some 

recursive techniques that have recently been developed by Ford and Sidi 
(ref. 6). These techniques are based on the determinant representations of 
s n k given in equations (2.18) to (2.20). For details we refer the reader to 
reference 6. 


As for the implementation of MMPE, this can be done by either solving the 
linear k x k system in equation (2.17) for the c ^ , i = 0,1, . . . , k-1, and 
then proceeding as in equations (2.4) and (2.5), or by using the recursive 
techniques of reference 6. 

The overhead in the implementation of MPE and RRE is largely due to the 
scalar products that are needed, c.f. equation (2.8) for MPE and equation 
(2.13) for RRE. To reduce this cost one might replace these inner products by 
a pseudo inner product involving only part of the components of the vectors 
x j . In computational fluid mechanics problems this could be done by excluding 
some of the dependent variables from the inner product. In discretized 
continuum problems, in general, one can also do this by excluding some of the 
mesh points from the inner product. 

The least expensive implementation in the above mentioned sense can be 
achieved by MMPE if the vectors qi in equation (2.17) are picked such that 
(qi, u n+j> require almost no computation. This can be achieved by taking the 
q^ to be the fundamental basis vectors, which amounts to picking k out of 
the N equations in equation (2.3) for determining c\, i =0,1,. . . , k— 1 . 
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One may, of course, consider other similar strategies, in which the vectors 
qi have very few nonzero components. 

Finally, we would like to mention some new extrapolation methods that have 
been proposed in reference 19 and are akin to MPE and RRE. These methods 
differ from MPE and RRE in the way the overdetermined systems of equa- 
tions (2.3) and (2.11) are solved. The standard linear least squares method 
minimizes the Q. 2 -norm of Uc + u n+ |<for MPE and of Wq + u n for RRE, where 
this norm is defined by ||z|| = y(z,z) with (y,z) as d efined in equa- 
tion (2.7). We can modify this norm to read ||z||m = y(z,Mz) for some 
hermitean positive definite matrix M. Going further, we T can choose to 
minimize the (weighted) ap-norms of Uc + u n+ |<_ and Wq + u n . In particular, 
the i]- and Q-co-norms give rise to problems that can be solved using linear 
programming techniques. 


3. CONVERGENCE AND STABILITY OF VECTOR EXTRAPOLATION METHODS 

We now state without proof some results concerning the convergence and 
stability properties of MPE, RRE, and MMPE. These results are helpful in 
understanding how these methods work, and how and when they can be used in an 
effective manner. All our results are stated for those sequences that arise 
from iterative solutions of linear systems of equations. 

Let s be the unique solution to the linear system 

x = Ax + b, (3.1) 

where A is a constant N x N matrix, and b is a constant vector in C^. 
Given xq, an initial approximation to s, generate the vectors x i , X£, 

. ... by the iterative method 

x j + .| = AXj + b, j = 0,1 (3.2) 

Let X], X£, . . . , X r be the distinct eigenvalues of the matrix A, and 
let V], V 2 , . . . , v r be corresponding eigenvectors. Obviously r < N. 
Assume for simplicity that A is diagonal izable so that the initial error 
xq - s can be expressed in the form 


x 


0 


- s 


5 vi 


(3.3) 


for some scalars a.. Consequently 

r 

x n - s = X a.v.X?, n - 1.2 (3.4) 

n i=i 1 1 1 

The assumption that equation (3.1) has a unique solution implies that \] * 1 
for all i. Without loss of generality we can further assume that (1) Xj * 0 
for all i, since a zero eigenvalue does not contribute to equation (3.4) for 
n > 1, and (2) ai * 0 for all i, since if <x\ = 0 for i = i ' , then we can 
discard it from equation (3.4) and rename the eigenvalues. 


I 


| 


! 
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Let us now order the \\ such that 

1*1 1 > M ^ 1^3 1 > • • • <3.5) 

We now state a convergence theorem for MPE, RRE, and MMPE, whose proof 
can be found in references 15 and 19. 


Theorem 3.1. Assume 


l*kl > l*k+ll- (3.6) 

Then, for both MPE and RRE, the approximations s n>k to s satisfy 

Sn , k. - s = r<n)|X k+ i|n, (3.7) 

where the vector T(n) is such that 

1 1 r<n) 1 1 =0(1) as n •* «, (3.8) 


and is proportional to n (X^ - 1) _1 . Furthermore, if we denote the y^ in 

i=l 

equations (2.5) by y^ n ’^ to stress their dependence on n and k, then 

k 

* - 1 1 \ 

K+ 1 1 as n co. ( 3 . 9 ) 


S >' ■ ri • i 


X, 


Equations (3.7) to (3.9) hold also for MMPE provided 

(q 1 ,v ] ) . . . (Q r v k ) 

(q k , v ] ) . . . (q,,^) 


k’ k 


* 0 . 


(3.10) 


Note that the results stated in Theorem 3.1 concern the strategy in which 
first a large number of iterations have been performed through equation (3.2) 
and then MPE or RRE or MMPE had been applied with fixed k. As such. Theorem 
3.1 has the following important implications: 

(1) If MPE or RRE or MMPE is applied to the vector sequence xq.xi, 

. . . , beginning with x n , for large n, then provided (3.6) holds, the error 
in s n k , the desired approximation, behaves like |X k+ ]| n . In view of the 
fact tfiat the error in x n+k+ i behaves like |X)| n for n large, and 
|Xi | > |x k+ ] | , all three methods accelerate the convergence of the sequence 
xo,xi,X 2 , .... when the latter converges. Furthermore, even when this 
sequence does not converge, i.e., | X-| | >1, s n k converges to s as n 
becomes large, provided |X k+ i| <1. 

In some applications the matrix A may have several distinct dominant 
eigenvalues that are all equal to the spectral radius p(A) in modulus. In 
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order to achieve acceleration in such cases we need to take k at least equal 
to the number of these eigenvalues. 

(2) If some of the are very close to 1 as is the case in many prob- 
lems of interest, then the quality of s , may deteriorate, as the error in 

k ' 

s . is proportional to n (X. - 1) for n large. In fact, the compu- 
n , k i =1 1 

tation of s . may become numerically unstable in this case as is seen from 

n f In. 

equation (3.9). Specifically, the coefficients of the polynomial 
k 

II [(X - X. )/(l - X.)] are large when some of X ] , . . . , X^ are very close 
^ = ^ ( n k ) 

to 1 in the complex plane. Thus the y. ’ , being approximations to these 

coefficients, are large too, although their sum equals 1, c.f. equation 
(2.6). This means that errors (round-off, in general) in the vectors xi are 
magnified severely in the computation of s n> ^, resulting in instability. One 
suitable way of overcoming this problem is to apply MPE, RRE, or MMPE to the 
sequence yj = x p j, j * 0,1, .... p being an integer greater than 1. With 
this equation (3.4) is now replaced by 

y n - s = £ a.v.p" , n * 1,2 (3.4)' 

n 1=1 1 1 1 

where p. = xP , i = 1 ,2 By picking p large enough we can cause p 1 , 

. . . , P k to be far from 1 in the complex plane, thus stabilizing the converg- 
ence acceleration methods. It is interesting to note that equation (3.2) 
becomes 



3 - 0 . 1 , 


(3.2)' 


(3) The already computed y j , 0 < j < k, can be used for estimating X], 
X 2 , . . . , X^, the largest eigenvalues of A. In fact, the zeros of the poly- 

k i 

nomial £ y.X J are the estimates of interest. If we denote these zeros by 
j=0 J 

X^n), . . . , X^(n), then with appropriate ordering 

IX, 


X.(n) - X. = 0 


k+1 


'j 


as n + 1 < j < k. 


(3.11) 


see reference 17. 


We note that the results of Theorem 3.1 and equation (3.11) do not hold as 
they are, in general, when the matrix A is not diagonal izable, although simi- 
lar but slightly complicated results for this case exist. For the precise 
statements and proofs of these results see reference 18. 
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As has already been mentioned, Theorem 3.1 above explains the convergence 
and stability properties of MPE, RRE, and MMPE when k is held fixed and n 
is increasing. Another obvious use of these methods is one in which n is 
kept fixed and k is increasing. An error analysis for this use, pertaining 
to MPE and RRE, has been given in reference 16. Part of this analysis is sum- 
marized in the following theorem. For simplicity we fix n = 0 and consider 
the sequence s^ = so,k. k = 0,1,2, . . . , with sq = sq q - *0- 


Theorem 3.2 Let the matrix C = I - A be such that Cf, = (C + C*)/2, the 
hermitian part of C, is positive definite. (This implies that all 

eigenvalues of C have positive real parts.) Then for some appropriate norms 
11 . 11 " 



where L is a constant independent of k, dependent on the norm and the 
extrapolation method, and is the minimum of ||Qk<C)||, the natural 
matrix 8.2 -norni °f Qk (C) > over a11 polynomials Qk<z> of degree at most k 
satisfying Qk<0) =1. Let now F(d,c,a) be the smallest ellipse that 
contains all the eigenvalues of C and is itself contained in the right half 
of the complex plane, such that its center is at d, its foci are at d±c, and 
its semimajor axis length is a. Then in (3.12) can be bounded as 

r k < I3k m n\ (3.13) 

where |3 is a constant independent of k, m is a fixed nonnegative integer, 
and n is given by 


n = exp (4>-Rew) < 1 

( 

with 4> = cosh -1 (a/ 1 c | ) and w = cosh -1 (d/c). 

(When the matrix C(or A) is diagonal izable m = 0.) Consequently, we can 
replace (3.12) by 



Tk m „ k 




Connections with Krylov Subspace Methods 


(3.15) 


It is furthermore shown in reference 16 that MPE and RRE are Krylov sub- 
space methods and that they are equivalent to the Arnoldi method and the method 
of generalized conjugate residuals respectively, when the latter are being used 
for solving the equations Cx = b with the matrix C as in Theorem 3.2. 

Recall that when C is hermitian the Arnoldi method reduces to the method of 
conjugate gradients and the method of generalized conjugate residuals reduces 
to the method of conjugate residuals. In this sense then MPE is a generaliza- 
tion of the Arnoldi method, which in turn, is a generalization of the method 
of conjugate gradients. Similarly, RRE is a generalization of the method of 
generalized conjugate residuals, which in turn, is a generalization of the 
method of conjugate residuals. For details and the relevant literature see 
reference 16. Although MPE and RRE can be applied in a very straightforward 
way to the vector sequence obtained by a fixed point iteration technique, for 
nonlinear as well as linear problems, the other methods require some kind of 
local linearization for nonlinear problems that may increase the cost of 
acceleration depending on how the linearization is performed. 
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From what has been said in the previous section it is obvious that we 
would not be interested in increasing k indefinitely as this would require an 
increasing amount of storage. From Theorem 3.2 however, it is easy to see that 
the following strategy, which has been designated "cycling" in reference 20 is 
useful in solving a system of the form x = F(x): 

Step 1. Pick s£ 0) □ x Q arbitrarily and set q = 1. 

Step 2. Generate x ] x k+1 by x.. + 1 = F(x^), j = 0,1, . . . , k. 

Step 3. Use MPE, RRE, or MMPE to compute s^ q) from x^.x^ .... x k+1 . 

Step 4. If s^ is a satisfactory approximation, stop; otherwise, replace 

x 0 by ’ and q by q+1 ’ and 90 bo Sbep 2 ' 

If F(x) is a linear function of x, then Theorem 3.2 can be used to prove 

that 



This implies that taking k sufficiently large we can cause LT^ < Tk m r^ < 1- 
This, by (3.16), guarantees the convergence of cycling. (3.16) for this case 
also suggests that cycling is a linearly converging process. 

In case \] \ are very close to 1 we can modify Steps 2 and 3 

to read 

Step 2'. Generate yo, yi yk+1 by Xj + i = F(x-j), j =0,1, 

. . . , p(k+l ) - 1 , and yj = x p j, j = 0,1, . . . , k+1 . 

Step 3'. Use MPE, RRE, or MMPE to compute s^ q) from yo » y 1 . - - - » yk+1 • 
This helps to stabilize the extrapolation process. 

4. APPLICATIONS 

We shall now mention some applications to computational fluid mechanics 
problems. In this respect the following remarks are important. 

(1) Practically no acceleration is achieved for problems that use explicit 
schemes. The addition of even a small amount of implicitness, such as the 
implicit residual averaging (ref. 8), makes the scheme quite suitable for 
acceleration. We note that the implicit residual averaging is a simple device 
that acts like a filter on the approximate solution produced by the explicit 
scheme, and can be added to the scheme without making any changes in it. 

(2) The vectors that are processed by the extrapolation methods must 
include the boundary values as a rule. This is especially so for problems in 
which the boundary values are time-dependent. Going back to the discussion of 
our introduction the discrete operator in equation (1.3) also contains the 
boundary values when these are time-dependent. Thus the boundary values 
influence the yj in the extrapolation process, and are influenced by the yj . 
(Time-independent boundary conditions, however, neither influence, nor are 
influenced by, the yj in the extrapolation methods.) 
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(3) As they are described in the previous sections, the extrapolation 
methods operate on the vectors xg, , , . . , in a global manner, i.e., the 
y's in equation (2.5) are the same for all components of the vectors x j . In 
some applications we may want to apply these methods to different portions of 
the vectors Xj separately, and this may produce more acceleration in some 
problems. For computational fluid mechanics this may amount to accelerating 
each of the dependent variables separately. 


Example 1. Ducted Flow Over an Axisymmetric Center Body 

The duct is cylindrical in shape and the center body is the hub of that 
found on the General Electric unducted fan. The flow is assumed to be 
inviscid, and thus the Euler equations are solved. In addition, the flow is 
assumed to be axisymmetric. 

This problem was solved by using two three-dimensional codes that employed 
the same mesh. 

The first code uses an explicit finite volume four-stage Runge Kutta 
scheme patterned after Jameson, Schmidt, and Turkel (ref. 9). The dependent 
variables for this code are pu, pv, pw, p, pe 0 , where u, v, and w are the 
axial, radial, and circumferential components of the velocity in cylindrical 
coordinates, p is the density, and e© is the total internal energy per unit 
volume. The scheme also contains artificial viscosity. It is supplemented by 
local time-stepping and by implicit residual averaging in the axial and radial 
directions. For a proper description of this scheme see reference 4. 

The second code uses an implicit finite volume flux-vector splitting 
scheme. The dependent variables for this scheme are pu, pv, pw, p, and pe 0 
as in the first code, except that now u, v, and w stand for the components 
of the velocity in Cartesian coordinates. For a proper description of the 
scheme used in this code see reference 1. 

The mesh used by both codes consists of 56 points in the axial direction 
and 7 points in the radial direction. The first code also requires a minimum 
of 5 points in the circumferential direction as dictated by the fourth order 
difference employed in representing the artificial viscosity. We chose to have 
7 equally spaced points in the circumferential direction their spacing being 
2ir/56 rads. If periodicity is imposed in the numerical solution of an axisym- 
metric problem, then there should be no variation in the circumferential direc- 
tion. This was achieved to machine accuracy in the iterative scheme used in 
the first code. The mesh employed was generated algebraically as described in 
reference 13, and it is shown in figure 1. 

Both codes were run at a Mach number of 0.60, where the critical Mach 
number for this problem is approximately 0.62. 

In figures 2(a) to 3(b) there are six curves numbered 1, 2, .... 6. 
Curve number 1 is for the Q^-norm of the error associated with all five 
dependent variables. Curves number 2, 3, . . . , 6, are for the 0,2-norms of 
the errors associated with p, pu, pv, pw, and pe 0 respectively. 

Figure 2(a) shows the results obtained by using the first code without 
extrapolation. As is seen the residual history gets into the linear regime 


13 


after 1000 iterations. Figure 2 shows the results obtained by the first code 
with extrapolation being applied only once to the sequence x n , x n+ p, x n+ 2 p, 

• • • > x n+(k+l)p with n = p = 10, and k = 20. The approximation 
obtained by the extrapolation method is then used as the new x 0 and new vectors 
are obtained from it by iteration. We see that with only one extrapolation the 
0 . 2 -norm of the error has dropped from about 10~ 4 -3 to about 10 _ 6 - 8 ( resulting in 
a gain of about 2.5 orders of magnitude. The jump in the & 2 -norm of the 
residual vector for pw, to circumferential momentum per unit volume, is a 
k 

result of Iy-sI being relatively large. This can be explained as follows: 
j =0 1 31 

Theoretically pw is supposed to be zero. Because we are using a three- 
dimensional code we can achieve at best machine zero for pw, and this is the 
case for the first code. Thus the residual vectors associated with pw are at 
best machine zero. By equation (2.5) then pw and hence its residual are of 
k 

£ 

j=o 

ure 2 (c) shows the results obtained from the first code with extrapolation in 
the cycling mode starting at the 100 th iteration, with p = 10 and k = 20 
again. We again observe a substantial amount of acceleration. However, the 
amount of acceleration now is smaller than that we observe in figure 2 (b). 

The reason for this may be that we are applying the extrapolation method much 
before the residual history reaches the linear regime. The jumps in the 
2 - 2 -norm of the residual vector associated with pw can be explained as above. 

Figure 3(a) shows the results obtained by using the second code without 
extrapolation. Figure 3(b) shows the results obtained by the second code with 
extrapolation in the cycling mode starting at the 10 th iteration with p = 1 
and k = 10. We see from this figure that with the second code (involving 
implicit flux-vector splitting) cycling, even when started very early in the 
iteration process, is very effective in this case. A reduction of approxi- 
mately 9 orders of magnitude is achieved without extrapolation in 1000 itera- 
tions, whereas with extrapolation this takes only 500 iterations. Also machine 
zero is reached in about 1400 iterations without extrapolation and in about 800 
with extrapolation. 

Before closing we also mention that we assessed the rates of convergence 
of both codes by estimating the largest eigenvalues of the Jacobian matrix 
i|»'(u). We recall that this is done by solving the polynomial equation 

k i 

]T y.X J = 0 , where the y's are those associated with s n . for n suffi- 
j =0 3 n,K 

ciently large, c.f. equations (3.9) and (3.11). The modulus of the largest 

k ■ 

zero of the polynomial Y-^ J is an estimate for pty'Cu)], the spectral 

j =0 3 

radius of <j>'(u). The information on the largest eigenvalues of h*'^) is then 
also used in deciding what values to take for p and k. 


Example 2. Flow Over a Single-Blade Row Unducted Fan 

The geometry consists of an eight-bladed unducted fan designed to operate 
at a free-stream Mach number of 0.80 and advance ratio 3.1145. The hub for 
this fan is the one that was considered in Example 1. 



the norm of the residual for 
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We again assume the flow to be inviscid, thus we are solving the Euler 
equations in this case too. 

The code used in this example is the first code that was used in 
Example 1. However, this time the implicit residual averaging is used in all 
three directions. 

The mesh used by the code consists of 88 points in the axial direction, 

23 points in the radial direction, and 11 points in the circumferential 
direction across one blade pitch. Ten points lie on each blade axially, and 
11 points, radially from the hub to the blade tip. A sting whose diameter is 
equal to the sting at the hub exit is affixed to the front of the nacelle. 

The mesh is generated algebraically as described in reference 13, and is shown 
in detail in reference 4. 

The six curves shown in figures 4(a) and 4(b) and numbered 1,2,..., 

6, again correspond respectively to the 2 , 2 -norms of the residuals associated 
with all five dependent variables, with p, pu, pv, pw, and pe 0 respectively. 

Figure 4(a) shows the residuals obtained without extrapolation. As is 
seen the residuals drop by about 2.5 orders of magnitude in 1000 iterations. 
Figure 4(b) shows the residuals obtained by employing RRE in the cycling mode 
starting at the 100th iteration with p = 10 and k = 20. The amount of 
acceleration achieved by doing this is quite substantial, in that in 1000 
iterations the residuals drop by 4.5 orders of magnitude. 
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FIG. 1 AXISYMMETRIC CENTER BODY GRID FOR EXAMPLE 1. 
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